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ABSTRACT 

In this work we study heat and mass transport, fluid motion and solid/liquid 
phase change in the process of steady Bridgman growth of Pb. 8 Sn 2 Te 
(LTT) in an axially-imposed uniform magnetic field under terrestrial and 
microgravity conditions. In particular, this research is concerned with the 
interrelationships among segregation, buoyancy-driven convection and 
magnetic damping in the LTT melt. The main objectives are to provide 
a quantitative understanding of the complex transport phenomena dur- 
ing solidification of the non-dilute binary of LTT, to provide estimates of 
the strength of magnetic field required to achieve the desired diffusion- 
dominated growth, and to assess the role of magnetic damping for space 
and earth based control of the buoyancy-induced convection. The prob- 
lem was solved by using FIDAPf and numerical results for both vertical 
and horizontal growth configurations with respect to the acceleration of 
gravity vector are presented. 


N 

f NASA does not endorse commercial products. Details about the products named in this paper were 
included for completeness and accuracy. No endorsement or criticism of these products by NASA 
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1. INTRODUCTION 

Direct, quantitative analysis of convection and solute segregation in molten metals 
and semiconductors is particularly difficult. These materials are generally opaque, which 
hinders non-intrusive measurements. Thus, most studies axe limited to indirect measure- 
ment techniques. The role of numerical modeling then becomes crucial in analyzing the 
system. A better understanding of the complex relationship among thermal conditions, 
natural convection, growth morphology, and macrosegregation in solidification process can 
be achieved through synergistic theoretical [1,2] and numerical analysis [3,4] in combina- 
tion with various experiments on earth and in space. In numerical modeling, the geometry, 
furnace temperature settings, materials and other parameters can be easily altered and the 
impact of the changes in design parameters can then be examined in a more efficient and 
economic way. The full-scale modeling combined with the scaling and asymptotic anal- 
ysis provides a convenient and useful tool for the furnace design and growth condition 
optimization. 

Directional solidification, such as Bridgman growth, has been widely used in space 
experiments and in many industrial processes for electronic materials. As one of the most 
fundamental arrangements, Bridgman growth is of significant technological importance, 
both from experimental and numerical points of view. Therefore, in this paper we shall 
focus on the modeling of Bridgman growth. However, once validated, the model and numer- 
ical techniques developed in this work are generally applicable to a number of solidification 
processes. 

Dining crystal growth, temperature and concentration gradients often induce natural 
convection. It has been shown both experimentally and numerically that the flow mode, 
the shape of the solid/liquid interface and the thermal and solute profiles in a solidifying 
liquid are greatly affected by natural convection. The variation of convective strength has 
a direct impact on solute distribution (segregation). The distortion of the interface, in 
turn, affects convection. Thus a complex interaction exists among convection, thermal 
and solutal gradients and interface shape. 

In the vertical Bridgman growth, the furnace is assumed to be aligned perfectly parallel 
to gravity vector. This arrangement results in a great simplification for numerical treat- 
ment. In this case, the flow, thermal and solutal fields may be taken to be axi-symmetric, 
which is valid for laminar flow and axi-symmetric boundary conditions. When convection 
gains enough strength, the flow on the axi-symmetric plane forms one or multiple torus-like 
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cells [5], depending on the local thermal and solutal gradients. At 1-g level, the magnitude 
of the convective flow may be very high, resulting in time-dependent or even turbulent 
flows requiring transient 3-D analysis [6,7]. 

Another fundamental Bridgman arrangement is the horizontal configuration in which 
the furnace axis is aligned perpendicular to the gravity vector. As shown by Arnold 
et al [8], the flow mode during 1-g horizontal growth is the so-called shallow cavity flow 
with a complex 3-D structure. During shuttle flight in space, however, the magnitude 
and orientation of the gravity vector are generally a priori unknown functions of time 
(unsteady). Consequently, the flow structure and segregation fields in the solidifying liquid 
during space growth are still theoretically unknown. 

In the solidification process, solute concentration at the interface is affected by the 
removal rate of excess solute rejected at the moving interface. Under ideal conditions (no 
convection), the rejected solute is transported by diffusion only. In reality, however, the 
transport of excess solute to and from the interface is greatly affected by additional factors 
such as the natural convection. In this study we consider buoyancy forces induced by both 
thermal expansion and composional differences. For PbSnTe, it is important to note that 
the magnitude of density variations due to thermal and solute gradients are globally of the 
same order and a positive change of AT and AC will reduce the density in the buoyancy 
term in the momentum equation. The segregation coefficient k<l and p S oiute < Psoivent 
imply that lighter solute is rejected at the interface. Consequently, the double-diffusive 
convection is solutally unstable in a classical vertical growth system. 

The conventional segregation quantity provided by experimental measurement is the 
effective segregation coefficient defined as the ratio of the concentration in the crystal to 
that in the bulk of the melt. Due to the existence of a solute boundary layer adjacent to 
the solid/liquid interface during the solidification process, this quantity differs from the 
interface segregation coefficient, which is calculated by the ratio of solute concentration in 
the grown crystal to that at the liquid side of the interface. The difference between these 
two quantities will depend on the structure and strength of the flow in the melt as well 
as on the growth rate. When flow is dominated by natural convection, the relationship 
between these quantities and the complex interactions among the driving forces can be 
obtained only from a full numerical simulation [9]. 

The idea of using a magnetic field to damp melt turbulence and thereby improve mi- 
croscopic homogeneity of the crystal was introduced in 1966 independently by Utech and 
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Flemings [10,11] and by Chedzey and Hurle [12]. It was shown in their work that a trans- 
verse magnetic field damped temperature fluctuations in the horizontal directional solidifi- 
cation of tin and InSb. In 1981 workers at Sony in Japan [13] developed a growth process 
using a static transverse magnetic field generated by an external electromagnet. Since then 
there has been considerable amount of research efforts pertaining to the magnetic damping 
in the growth of silicon, the H-V compound semiconductor, gallium arsenide, and indium 
phosphide, etc. [14]. Motivation for the use of an applied magnetic field has been widened 
subsequently. In addition to damping out the chaotic and time-periodic convection in the 
melt and consequently reducing fluctuations in solute concentration, magnetic field can 
also be used to control the growth conditions at various stages in the growth process. For 
instance, the alteration of the laterally averaged axial concentration of oxygen is one of the 
primary advantages of using magnetic damping in Czochralski growth of silicon [15,16]. 

There are several different configurations of the applied magnetic field reported in the 
literature. The axial and transverse fields are perhaps the simplest configurations am ong 
them. In addition, the concept of a configured field [17,18] was also developed to tailor 
the field configuration to the flows in the melt so that only the harmful flows are damped 
while the beneficial flow patterns are still retained. One such configured field is the cusped 
field formed with a pair of Helmholtz coils operated in opposed-current mode. The most 
extensive work to date, however, has been on the axial magnetic field. 

The application of steady-state magnetic fields offers a practical means of suppressing 
natural convection due to both the steady and unsteady changes in the residual accelera- 
tion. When an axial magnetic field is imposed, its effect on the convection in the melt is to 
interfere with the radial velocity component. As shown by asymptotic analysis [46,19], the 
magnitude of the radial velocity decreases proportionally to the square of B, the strength 
of the magnetic field. When B is sufficiently intense, an almost uniaxial flow, and hence 
the desired diffusion-dominant growth condition, may be obtained. 

In the literature, previous studies on vertical Bridgman growth with magnetic damping 
include both experimental and numerical efforts. K.M. Kim [20] grew InSb with a magnetic 
field greater than 1.7 Kilo- Gauss (KG) and obtained a striation free boule. The most 
obvious effect of an increasing magnetic field was the reduction and then the elimination of 
temperature fluctuations in the melt. Sen et al [21] found that a 2 KG transverse magn etic 
field was sufficient to statistically reduce the number of area defects in InGaSb, but was 
insufficient to produce an observable effect on the composition distribution. Matthiesen 
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et al. [22] demonstrated melt stabilization in doped germanium when grown in a 30 KG 
axial magnetic field. They observed a transition from convection dominated to diffusion 
controlled growth, for the dopant atoms, via the application of the magnetic field. Su, 
Lehoczky and Szofran [23] grew HgCdTe and HgZnTe in a 5 KG transverse field. Their 
results show a strong coupling between convection and the magnetic field but they were 
unable to fully suppress convection. 

In the pioneering numerical magnetic interaction studies for vertical Bridgman growth 
[19], D.H. Kim et al performed a finite element analysis for the growth of dilute Ga doped 
Ge. Their study shows that the compositional radial non-uniformity is greatest for an 
intermediate field strength. Stronger fields suppress convective flow completely and lead 
to uniform solute segregation across the crystal and to diffusion-controlled axial segrega- 
tion. They calculated that the effect of the magnetic field was very similar to reducing 
gravity for vertical growth. Motakef [24] has analyzed the growth of Ge, GaAs and CdTe 
in both low gravity and magnetic fields. He computed the ranges of ampoule sizes and 
gravity level or field strength where these single component materials could be grown under 
diffusion control. In a more recent work, Prescott and Incropera [25] simulated magnet- 
ically damped convection in solidifying Pb.igSn alloy. Their results show that magnetic 
damping significantly affects thermally driven flows during early stage of solidification. 
However, interdendritic flows and macrosegregation patterns are not significantly altered 
by moderate magnetic fields. Their scaling analysis suggests that extremely strong fields 
would be required at 1-g to effectively dampen convection patterns that contribute to 
macrosegregation. Simulations in their work were based on a continuum model for den- 
dritic solidification systems [26]. In [27], their work was extended to account for turbulent 
flow and to application of a time varying magnetic field which augments thermal buoy- 
ancy forces in the melt while opposing solutal buoyancy forces in the mushy zone. Due to 
the complexity of the problem and the large number of parameters involved, it is difficult 
to obtain the whole picture of the physical phenomena from numerical simulation alone. 
Hjellming and Walker [47-49] have therefore used asymptotic analysis to deduce the nature 
of the boundary layers formed. 

Among the previous efforts cited above, only two studies [21,23] involve the melt stabi- 
lization of a non-dilute multicomponent system. One important consideration for analysis 
of non-dilute alloys is the need to include the density variation due to compositional dif- 
ferences. For non-dilute alloys, our results (and the results in [5]) show that the effect of 
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the compositional derived density differences is important and should not be neglected. 
Furthermore, the non-trivial interaction between the thermal and the solutal fields and 
their coupled effects on the flow field in the melt through the buoyancy force term intri- 
cately depend on the magnitude and orientation of the gravity vector. Therefore we c ann ot 
decipher the interactive thermal and solutal effects a priori by order-of-magnitude analysis 
of the non-dimensional groups alone. 

In this study, we extend previous works in this area to non-dilute multicomponent 
systems and to the combined effect of magneto-hydrodynamic (MHD) and low gravity 
stabilization. In particular, this research is concerned with the study of the effects of an 
axial magnetic field on the Bridgman growth of Pb.gSn. 2 Te with convective flows induced 
by both thermal and solutal buoyancy forces. To this purpose, a complete 2-D (including 
axi-symmetric) finite element model (for full-coupled flow, thermal and solutal fields) has 
been established. The model considers the coupled momentum, energy, and mass transport 
equations with the Boussinesq approximation applied to the temperature and concentration 
buoyancy terms. The influence of the magnetic field on the flow is expressed through 
the Lorentz force. The temperature of the solid-melt interface is determined from the 
phase diagram. The interface position is simultaneously solved using the front tra cking 
approach in conjunction with a deformable mesh. The furnace temperature profile is 
imposed through a heat transfer boundary condition on the outer surface of the cartridge, 
details of which were determined from previous experimental data. The so-called pseudo- 
steady-state-model [5,19] is used to simulate the steady growth in both axial and horizontal 
Bridgman configurations. The problem is solved by the co mm ercial code FIDAP [28]. 

The primary objectives of this study are to provide a quantitative underst anding of 
the complex transport phenomena during solidification of non-dilute binarys, to furnish a 
numerical tool for furnace design and growth condition optimization, to provide estimates 
of the required magnetic field strength for low gravity growth, and to assess the role of 
magnetic damping for space and earth control of the double-diffusive convection. These ob- 
jectives will be achieved via a systematic examination of the heat and mass transport and 
fluid flow phenomena using both earth and reduced gravity conditions, and by using MHD 
damping at various gravity levels. As an integral part of a NASA research program, our 
numerical simulation collaborates and supports both the flight and ground-based experi- 
ments in an effort to bring together a complete picture of the complex physical phenomena 
involved in the crystal growth process. 
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In order to provide a solid theoretical foundation for this work, this paper will be 
devoted to the modeling studies, leaving a more detailed description of the experimental 
work for subsequent reports. 

2. EXPERIMENTAL BACKGROUND 

The experiment [29,30] modeled in this paper studies the effect of the gravitational 
body force on the convective properties of the alloy compound semiconductor, PbSnTe 
(LTT), with body forces modified by both reduced gravity and by MHD. PbSnTe is an 
ideal material for this study because it was used in both a past flight experiment and 
a planned AADSF (Automated Advanced Directional Solidification Furnace) experiment. 
Both of these experiments are without magnetic field. Subsequent experiments using MHD, 
including earth based and in space, will form a comprehensive set of space processing 
experiments which will help to elucidate the gravity dependent physical phenomena for 
the growth of this class of materials. 

The bulk growth of this material is interesting from a purely scientific point of view, 
because the liquid is always solutally or thermally unstable. For the vertical growth, one 
of the pertinent fields, either temperature or concentration, will be in a stable orientation 
and the other field will be in an unstable orientation. This double convective instability 
cannot be made stable by simply balancing thermal and solutal expansion with a high 
temperature gradient [31]. 

PbSnTe is a semiconductor material with a compositional dependent energy bandgap 
which is adjustable from 0 ev («40% SnTe) to 6.4 ev (« 100% PbTe). It is a direct bandgap 
material hence it can be used for both infrared lasers and detectors. The utilization of this 
and other materials in this general class is dependent on both the crystalline perfection of 
the material and the compositional homogeneity. 

The first flight experiment, without magnetic field, in this series was on STS 61A 
in October 1985 [32]. It was shown that a high degree of convection was still apparent 
at the low growth rates used. The convective instability may have been exacerbated 
because the crystal axis and the steady state residual acceleration vectors were (estimated) 
perpendicular. The analysis of this flight sample served as a basis of continued ground 
based research and the development of the existing flight experiment. Previous studies 
show that the growth of PbSnTe is totally dominated by gravity driven convection on 
earth and is still strongly influenced by convection in space growth. These two facts make 
PbSnTe an ideal material for studying the effects of MHD damping on earth and in space. 
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The experiment consists of two parallel programs, the ground-based growth in the 
superconductor magnetic furnace and the space growth in a magnetically damped furnace. 
The ground-based experiments are conducted at MSFC where a crystal growth furnace 
with a 50 KG axial magnet and another furnace with a 5 KG transverse magnet have been 
developed. After growth, the samples are evaluated for compositional uniformity and defect 
structure. The compositional profile is the most sensitive measure of convection intensity 
during growth and is measured by a wavelength dispersive electron micro-probe. The 
planned flight experiment has basically the same design as that of the AADSF experiment 
but with the addition of the stabilizing magnetic field. 

The typical experimental parameters, such as the ampoule dimensions, temperature 
gradient, etc., are listed in Table 1. The ampoule materials are quartz and Inconel™. 
The thermophysical properties of PbSnTe can be found in refs. [29,33-36]. The phase 
Diagram is given in [37]. The electrical conductivity can be found in [38]. Other material 
properties are given in [39]. A typical thermal profile measured along the center of the 
ampoule is shown in Fig. 1. This profile is used as the reference temperature in the heat 
transfer boundary condition imposed in the mode ling 

3. THEORETICAL BACKGROUND 
3.1 Governing Equations 

In this paper, we consider heat and mass transport, fluid motion and solid/liquid phase 
changes in the crystal growth process. In particular, the liquid pseudo-binary mixture of 
LTT is assumed to behave as a Newtonian fluid and its motion is described by the following 
Navier-Stokes equation: 

to (^+“- Vu ) = - Vp + V • {p(T)[Vu + (Vuf)]} (la) 

+ Pog[l — fit(T — T 0 ) — fi c (C - Go)] + a m (u x B) x B 

where u is the velocity vector, p 0 is fluid density at the reference temperature T 0 , p is 
pressure, p is viscosity, T is the temperature variable, g is the acceleration of gravity, fit 
and fi c are thermal and solutal volumetric expansion coefficients. The Boussinesq model is 
adopted to approximate the buoyancy force caused by density variation with both temper- 
ature and concentration. A uniform magnetic field, B, is imposed in the axial direction of 
the furnace. When the ampoule and crystal are moved parallel to B, the imposed magnetic 
field does not cause a current in the materials. The electric field is everywhere equal to 
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Table 1. Experimental Design Parameters for the Vertical Bridgman Growth 


Parameter 

Symbol 

Value 

Ampoule Length 

L 

11 cm 

Sample Radius 

R 

0.5 cm 

Ampoule Outer Radius 

R a 

0.7 cm 

Thermal Gradient 

dT/dz 

«80°C7/cm 

Growth Velocity 

Vg 

lcm/hour 

Hot Zone Temperature 

Th 

1150 ° C 

Cold Zone Temperature 

T c 

550° C 

Temperature Difference 

> 

II 

£ 

1 

600K 



1200 

1100 

1000 

900 

800 

700 

600 



* 


* 


* 


* Tweak-3 Data 


* 


* 


* 








3K3K 




500 


J 1 1 I i L 


- 10.0 


- 5.0 0.0 5.0 

Z (cm) 


10.0 


Figure 1. Temperature profile obtained using a quartz sample (Cal9). It was mea- 
sured along the sample centerline with the furnace set to 525/575/1150/1150/1150 
and the sample translated at a speed of 2 cm/hour. 
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zero and the magnetic induction is not distorted by convection in the melt. Therefore the 
only effect of the magnetic field is the Lorentz body force term in eqn. (1), where <j m is 
the electrical conductivity of the melt. 

The LTT liquid is assumed to be incompressible and the incompressibility condition 
is given by 

V • u = 0 . (2) 

The heat transport is controlled by the balance of thermal energy 

Poc P + u • VTj = V • [«(T)VT] (3a) 

where c p is specific heat and k is thermal conductivity. The solute transport is governed 
by the balance of species concentration 

^ + u . VC = V • (DVC) (4a) 

where C is the concentration of impurity and D denotes the mass diffusion coefficient. On 
each segment of the boundary, it is necessary to prescribe appropriate boundary conditions. 
Details of the computational boundary conditions used in the modeling will be given in 
section 4. 

The change of phase from liquid to solid is described mathematically by the following 
phase conditions 


T l (S,t)=T s (S,t)=:T m , 

(5a) 

kiVTi ■ n — k s VT s • n = p s AH(u s — u/) • n , 

(6a) 

pi(\H - uj) • n = p s (u s - uj) • n , 

(7) 

(ui - u s ) x n = 0 

(8) 


which need to be satisfied at the solid/liquid interface. Here subscripts, s and l, refer to 
the solid and liquid region, respectively; u s is the solid pulling velocity; uj is the velocity 
of the interface; T m is the melting temperature; n is the unit norm of the interface pointing 
from the liquid to the solid; A H is the latent heat of fusion. From physical point of view, 
equations (5)-(8) represent thermal equilibrium at the interface, the heat flux balance 
between the phases which includes latent heat release, the mass flux balance across the 
interface and the no-slip condition at the liquid side of the interface, respectively. 
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In order to include the solute (species) transport in the phase change analysis it is 
necessary to add two more interface conditions: 

DtVCi • n - D S VC S ■h = C s {u s -u I )-h- C t (ui - uj) • n , (9) 

C s = kC t (10) 

where k is the partition coefficient. Eqn. (9) describes the mass conservation for solute 
transport across the interface and eqn. (10) is actually the chemical equilibrium deter- 
mined by the phase diagram. In this case the melting temperature depends on solute 
concentration and condition (5a) becomes 

T t = T s = T m (C) = T a + mCi (5b) 


where Ta is a constant and 

m = m(C) = ( n ) 

is the rate of change of the melting temperature with respect to C. Note that the model 
summarized in this section assumes a sharp solid/liquid interface and the problem will be 
solved by the interface tracking approach with deformable grids. 

3.2 Non-Dimensionalizadon 

It is well known that non-dimensional formulation renders many advantages for study- 
ing the physical problems. Scaling the fundamental variables with respect to their char- 
acteristic values and defining dimensionless parameter groups provide a measure of the 
relative importance of the various terms in the governing equations and identifies the dom- 
inant physical phenomena. For nonlinear problems, non-dimensionalization can also help 
ease the convergence difficulties in the nonlinear iterative solution process. For these rea- 
sons, the non-dimensional formulation is also used in this work, especially for the cases 
when the nonlinearity is high as well as for the parametric and asymptotic analysis. 

To non-dimensionalize the momentum equation, we first choose a characteristic veloc- 
ity U and a characteristic length L. Then we define the following dimensionless variables 


u* = u/U, x* = x/L 

p- = p/(poU 2 ), t* = tU/L 

T* = (r-T 0 )/AT, ^(T*) = /i(T)/#io 

Pi = 1 , Pt = PJ(PtAT) 


( 12 ) 
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Substituting (12) into (la) and introducing the non-dimensional parameters summarized 
in Table 2 leads to the following non-dimensional equation of motion 

/du* \ 1 

(^+u*-V*u*J = - VV + =V* • {fx*(T*)[V*u* + (V*u*) T ]}+ (16) 

1 HaP 

+ - <*“" + W*’ + x eB) x 6B 

in which e g and e b are unit vectors of gravity and magnetic field. 


Dimensionless Number 

Definition 

Range of Value 

Reynolds No. 

Re = pqUL/hq 

8.5 ~ 850 

Grashof No. 

Gr = /%/3rATgL 3 /fi% 

72 ~ 72 x 10 4 

Hartman No. 

Ha = BLy/<j m /nQ 

258 (5T) 

Thermal Peclet No. 

Per = UL/a o 

0.91 ~ 91 

Prandtl No. 

Pr = fj,QCpf kq 

0.107 

Thermal Rayleigh No. 

Ra T = pogfrATtf/noao 

7.7~7.7xl0 4 

Stefan No. 

St = AH/cpAT 

3.327 

Schmidt No. 

Sc = Pq/pqocc 

31.56 

Solute Peclet No. 

Pe c = UL/a c 

2.7~2.7xl0 4 

Solute Rayleigh No. 

Ra c = pog/3 c C 0 L 3 / p 0 a 0 

12 ~ 12 x 10 4 


Table 2 Non-Dimensional parameter groups and their typical ranges of value for vertical Bridge 
growth of LTT under the terrestrial and microgravity conditions. The upper limit is for lg and 
the lower bound corresponds to 10~ 4 g. T he characteristic velocity used here is based on the 
natural convection, i.e. U = TgL, where AT = Th— T m ~ 250° C. The characteristic 

length is the radius of ampoule, R = 0.5cm, which is different from the definition used in [19]. 

If the ampoule length, L = 11cm is used then some of the numbers will have much higher 
values. For instance, Rax will become as high as 8.2 X 10 8 . 

For heat transport, we have the following non-dimensional energy equation 

^ + u* • V*T* = * [a*(T*)V*T*] (3c) 

dt* Pe t 


Here the thermal diffusivity a*(T*) is assumed be a function of temperature and is scaled 
by a reference value qq- For constant thermal diffusivity, a* = 1. The dimensionless 
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convective heat transfer boundary condition reads 

Q* = h*[T* - T*) 

where 

n ._ «* ar . r. - ro 

W “ Pe T an*’ ' MV’ ° AT 

and the Stefan condition becomes 


(13a) 

(136) 




Qt =AH*(u:-u*j)-h 


( 66 ) 


where the dimensionless latent heat is given by 

AiT* = 1 


(6c) 


/t?o Ste 

Similarly, we have the following non-dimensional balance equation for concentration 


fi(~> i 

_ +u *.v*c=— v*.(a;v*o 


(4b) 


Note that C is already dimensionless and a* = a c /a° is the non-dimensional solute dif- 
fusivity. If a c = a Q c is a constant, then a* reduces to 1. The mass transfer boundary 
condition is transformed to 


q*m = h m C 


(14a) 


where 


— 

9m — 


a 


dc 


h* m = 


Pe c dn* 

and the phase change condition (liquidus) becomes 


paoPe c 


(146) 


T* = T^ + m*C (56) 

where = (T m - T 0 )/AT and m* = m/AT. 

3.3 Pseudo-Steady-State Model 

To simulate the steady growth of LTT in the vertical and horizontal Bridgman con- 
figurations, the so-called pseudo-steady-state model (PSSM) [5,19] is adopted in the present 
work. In PSSM, the directional steady movement of the solid/liquid interface during the 
steady growth is modeled by letting melt enter at its hot end with a uniform growth ve- 
locity V g and composition Cq and by removing the crystal from the cold end at a speed 
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that conserves the mass of the alloy in the system. Note that as a simplification of the 
real growth condition the PSSM neglects the transient effects in the field variables (such 
as velocity, pressure, temperature and concentration, etc.) caused by the steady decrease 
of the length of the melt during crystal growth and the displacement of the ampoule in 
the furnace. This simplification is valid for long ampoules and melts with low Prandtl 
numbers, in which the transient effects on heat transfer are small [40]. It is also known 
that for sufficiently long ampoules the thermal end effects are negligible [41]. 


Interface 


ir' 

J 

/ 

z'.z , 

z' 

z 

o , o 

► ’ o' 

►- o 

►- 

Liquid 

Solid 




h V* *h-*- Vg 

(a) time = 0 (b) time > 0 


• P{z,r) 


Figure 2 Definition of the moving coordinate system zor used in the PSSM and its 
relation with a reference frame z'o'r' fixed in space. 


Perhaps the most important feature of PSSM is the application of the moving coordi- 
nate system. In this work, the moving coordinate system is introduced by fixing its origin 
at the center of the moving solid/liquid interface. Fig. 2 shows the definition of the moving 
coordinate system and its relation with a fixed frame in space. Let zor denote the moving 
coordinate system for the axi-symmetric case and z'o'r' a fixed coordinate system in space, 
respectively. We choose the origins of the two systems to coincide at time equals zero, as 
depicted in Fig. 2(a). At £>0, the coordinates of point P in the moving coordinate system 
are (r, z) while its coordinates with respect to a stationary observer sitting on the fixed 
frame are (r', z'). From Fig. 2(b), it is evident that 


r = r 


z = z'~ Vgt 


( 12 ) 


Here V g is the growth rate with V g >Q corresponding to melting and V g < 0 to solidification. 

The relation of velocity in the two systems can be easily derived by considering the 
time derivatives of (12), namely 


r = r' u r = u' r 


(13) 
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z = z' -V g ' =*► u z + V g = u' z (14) 

According to (13) the r-components of velocity in the two systems are identical. Since the 
s/1 interface becomes stationary in the moving coordinate system, the steady movement 
of the interface is now represented by the translation of the moving coordinate system in 
an opposite direction to the growth (or melt) direction. The translation of the moving 
coordinate system does affect the ^-component of velocity. Equation (14) indicates that 
although the material particle in the solid is stationary with respect to the fixed coordinate 
system, it has a constant translation velocity when viewed from the moving coordinate 
system. Consequently, the velocity boundary conditions for u z need to be modified in the 
moving coordinate system. Let Cp denote the given velocity on boundary in the fixed 
coordinate system, viz. 

u' z - <p(r',z r ) . (15) 

Then the velocity boundary condition in the moving coordinate system becomes 


it, = z) - V g . V (r, z) £ T„ 


(16) 


For instance, at the no-slip boundary (p = 0, but under system zor eqn. (16) gives u z — —V g , 
which is due to the translation of the moving coordinate system. 

In order to obtain the transformed governing equations in the moving coordinate sys- 
tem, we also need to consider the spatial derivatives. By using the chain-rule, it is easy to 
verify the following relation for the first-order spatial derivatives 


d 

d dr 

__ d_ 

dr’ 

dr dr’ 

dr 

d 

d dz 

__ d 

dz ’ 

dz dz’ 

dz 


(17) 


And likewise for the second-order derivatives we have 

d 2 _ 

dr 12 dr 2 ’ dz 12 dz 2 

Eqns. (17) and (18) suggest that the spatial derivatives in the governing equations are 
unaffected under the coordinate transformation. Applying (12)-(14) and (18) in eqn. (3), 
the transient energy equation is now reduced to the following steady state form under the 
moving coordinate system 



p 0 c p (u + V)-VT = V-(kVT) . 


(36) 
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Here V is the translation velocity vector of the moving coordinate system. 

3.4 Special Treatment for Concentration 

The special difficulty in solving mass transfer in phase-change problems is the disconti- 
nuity of solute concentration itself across the interface, as indicated by interface condition 
(10). Furthermore, the melting temperature is no longer a constant, instead it has to be 
determined by the local solute concentration. However in the standard FEM formulation, 
both the trial and test function have to be at least C° continuous and no jumpls allowed for 
the primitive variables. Therefore a special treatment is necessary for the species equation. 
To eliminate the discontinuity of C at the interface, the following transformation 

c; = c./k ( 19 ) 

is performed for concentration in the solid region. The interface condition (10) is then 
reduced to 

Ci = c; ( 20 ) 

and thus C becomes continuous across the phase change interface. It can be proved [30] 
that the mass flux balance condition (9) can be rewritten in the following form 

{DiVCt - D S VC S ) • n = h c Ci (21) 

where the coefficient 

h c = p,(l - k)(n s - uj) • n + D s ( 1 - k)h • VC*fC t . (22) 

The second term in (22) involves spatial gradient of concentration at the solid side. Since 
the solute diffusivity in solids, i.e. D s , is usually several orders of magnitude smaller than 
Di, it is reasonable to assume this term is negligible. In our computation, the interface 
condition (21) is imposed through a species transfer boundary element at the moving 
interface. 

3.5 Dimensional Units for the Lorentz Force Term 

When the dimensional form of the momentum equation is used in the magnetic damp- 
ing analysis, it is important to have the units of the Lorentz body force term to be dimen- 
sionally consistent with the other terms in the same equation. It is easy to verify that the 
other terms in equation (la) have the following general units 

du force 
/3 ° dt volume 


(23) 
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For the standard international (SI) unit system, all the units involved in the Lorentz 
force term are clearly defined, such as the following definitions for the magnetic field 
strength and the electrical conductivity 

N 


T (Tesla) = 
1 


A ■ m 
A 2 • s 


(24) 

a-m N ■ m 2 ^ ^ 

where A is electic current in amperes, Q, is resistance in ohms and N is force in Newtons. 
Note that the following relations: O = V/A, V - W/A and W — N ■ m/s are used in 
deriving (25). By using (24) and (25) we can easily confirm the SI units of the Lorentz 
force term, i.e. 

A 2 • s 


'rrx 


b 2 u 


N • m 2 


f N V m _ 
\A • m ) s 


N 




(26) 


which is indeed consistent with general units (23) of the other terms. 

However, for the centimeter-gram-second (CGS) unit system, things are not so 
straightforward. This is because there are no CGS units commonly used for voltage, 
current, etc. In this case, we need to consider the following term 

(27) 


—B 2 = a m B 2 
Po 


where a m = <T m /p 0 is called the mass diffusivity in FEDAP. Then we introduce the following 
result 

Proposition: The dimensional unit of {27) is invariant for both the SI and the CGS unit systems. 
Proof: Substituting the SI units for cr m , po and B into (27) gives 


^rn jg2 


A 2 


s m 


3 


/ N V _ s • N _ 1 
\A-mJ kg-m s 


(28) 


po N • m 2 kg 
which is the same in the two unit systems. 

To confirm the above results, we substitute (28) into the Lorentz force term and obtain 

g 1 cm dyne 


Po(x m B U ~ — o = 


cm° s s 


cm c 


(29) 


which is exactly the consistent unit for the CGS unit system. Therefore, the above results 
suggest that for magnetic damping analysis based on the CGS units, one can use a mixed 
units in the Lorentz force term, namely use SI units for the magnetic field related term 
(27) while inputing the rest of the quantities in the Lorentz force term with CGS units. 
The CGS units of the Lorentz force term will then be consistent with the other terms in 
the momentum equation. 
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4. NUMERICAL SIMULATION 
4.1 The FEM Model 

In this paper we consider two Bridgman growth configurations. The first is the bot- 
tom seeded (vertical) Bridgman growth. In this configuration the hot (melt) zone is on 
the top and the cold (crystal) zone is at the bottom. The furnace axis is parallel to the 
gravity vector which is in the vertical direction pointing downwards. For axisymmetric 
boundary conditions it is reasonable to assume that the heat, species and flow fields are all 
axisymmetric, i.e. independent of the ^-coordinate. Therefore a simplified axisymmetric 
model, as depicted in Fig. 3, is used to model the vertical Bridgman growth. The compu- 
tational boundary conditions imposed for the energy, momentum and species equations are 
summarized in Figure 4. For the momentum equation, we specify the no-slip conditions 
on the surface between the sample and ampoule wall. For heat transport we impose the 
measured thermal conditions on the outer surface of the ampoule. As shown in Fig. 4(b), 
the thermal profile on the ampoule surface consists of the hot, adiabatic and cold zones. 
The concentration boundary conditions are shown in Fig. 4(c). 


■« Growth Rate V g ► gravity 


R 

z 

o *■ 

H 6.5cm + 4.5cm H 

Figure 3. Schematic diagram and geometric definitions of a simplified axisymmetric 
FEM model for the vertical Bridgman growth configuration. 

The second configuration is the horizontal Bridgman in which the gravity vector is 
perpendicular to the furnace axis (i.e. the 2 -axis defined in Fig. 3). In this case the 
solutions are assumed to be symmetric about the vertical center plane and an idealized 
two-dimensional model is used to simulate this center plane. Although this 2-D model is 
just a simplification for the real situation, we hope it can provide at least some qualitative 
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Adiabatic Zone 



dT/dR = 0 

(a) Thermal Boundary Conditions 



u T — 0 u z = free o Z 

(b) Velocity Boundary Conditions for Pseudo-Steady State Analysis 



dC/dR = 0 

(c) Boundary Conditions for Solute Concentration 


• Figure 4. Computational boundary conditions used for the axisymmetric model. 







20 


M. Yao, A. Chait, A.L. Pripp and W.J. Debnam 


information for this growth configuration. More accurate modeling will rely on the full 
three-dimensional model which is currently under development. 

4.2 The Front Tracking Approach 

The most challenging difficulty posed by phase change problems is the moving solid- 
liquid interface whose position is usually an unknown function of time and space and needs 
to be determined as a part of the solution. In the literature various numerical techniques 
have been proposed to deal with the moving interface for phase change problems [42], 
Among them, two main classes of methods can be distinguished, namely the fixed-grid 
enthalpy methods and the front-tracking methods. The front-tracking technique with 
deforming mesh is used in our modeling. 

Unlike the enthalpy method, the phase change front tracking method can model the 
phase change problems with a sharp (single) solid/liquid interface. The front tracking 
approach used in this work involves a deformable spatial mesh in which nodes located on 
the interface are allowed to move such that they remain on the moving boundary. For each 
node on the moving interface, an additional degree of freedom is introduced. This new 
degree of freedom determines directly the position of the node in space and is an integral 
part of the representation of the moving interface. 

To update the interface position and remesh the interior domains, a method of spines 
is used. It is a generalization of the method developed by Saito and Scriven [43], in which 
the moving nodes lie on and the interface movement is guided by the generator lines called 
spines. In particular, the position of the moving node is represented parametrically by 

Xi = a x [h j+ 1 + uti(hj - h J+ i)] + (3 X 

Vi = a y [hj+i + - hj+i)] + f$ y (30) 

Zi = a z [hj+i + u)ti(hj — hj+i)] 4- (3 Z 

where hj is the interface location parameter for a given spine, (a x ,a y ,a z ) is the direction 
vector, (fl x ,P y ,(3 z ) is the base point of the spine. The location (xj,yi,Zi) of the moving 
node on the spine is determined from its relative position, uti, to the moving interfaces 
located at hj and hj+ j . Here hj are the new degree of freedom introduced. 

For steady state problems, applying Galerkin’s formulation and the standard dis- 
cretization procedure to the momentum equation (1) results in a nonlinear algebraic system 
in the following matrix form [28] 


A(U)U + K(T, U) - CP + BX = F 


(31) 
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where X is the global vector of the moving interface unknowns, A(U) accounts for contri- 
bution from the convective terms, K(U) includes the diffusive terms, C is the divergence 
matrix, B represents the contribution of the normal stress balance condition at bo un dary, 
F is a body force vector. 

4.3 The Segregated Solver 

System (31) along with the energy equation, species equation and other interface and 
boundary conditions are solved by a segregated solution procedure. In contrast to the tra- 
ditional fully coupled approach (such as the Newton-type solver), the segregated solution 
algorithm does not form the global system matrix directly. Instead, it decomposes the 
global system matrix into smaller sub-matrices each governing the nodal unknowns associ- 
ated with only one conservation equation. These segregated sub-matrices are then solved 
in a sequential manner. Since the storage of the individual sub-matrices is considerably 
less than that needed to store the global system matrix, the storage requirements of the 
segregated approach are substantially less than those of the fully coupled approach. 

4.4 Numerical Solution 

The axisymmetric and 2-D FEM models are built with the 4-node bilinear element, 
in which velocity, temperature and species are approximated by bilinear shape functions. 
The pressure is approximated as piecewise constant. 

In this work numerical solutions are obtained by using a modified version of the finite 
element program FIDAP. The details of the FEM formulation in FIDAP are documented 
in [28]. The nonlinear iteration termination is controlled by a specified tolerance of 10 -3 
for the relative error norms of velocity, residual and free surface update. 

5. BENCHMARK TESTS 

Due to the complexity of the crystal growth problem considered herewith, there ex- 
ists no closed-form (analytical) solution, nor adequate published results in the literature. 
Therefore, it is necessary to conduct a series of benchmark tests to check the validity 
of the numerical model and the accuracy of the numerical solutions. In this section we 
summarize some of the typical tests that have been performed. These test problems are 
usually simplified so that analytical solutions are available for quantitative comparisons. 
Each simple test is aimed to examine only a particular part of the model. 
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5.1 1-D Crystallization of A Binary AUoy 

First, we consider the classical 1-D solidification of a binary alloy in a semi-infinitive 
region, r >0. Initially the whole block of molten alloy is in liquid phase with uniform initial 
temperature and concentration distributions, To and C 0 . At t = 0 + , the temperature at 
the left end surface, x = 0, is lowered to T\ < T m . Hence crystallization of the mixture is 
initiated at the boundary a?=0. The complete mathematical formulation for this problem 
can be found in [44]. The main objective of this test is to verify the formulation and FEM 
implementation related with concentration, phase change conditions and front tracking. 

Under certain assumptions, such as constant conduction and diffusion coefficients 
in each phase, this problem admits a similarity solution which is usually referred to as 
Rubinstein solution and has the following general form 


S = 2pVt 

C s = k s T m 

T.=Ti+ (T m - Ti)erf M (£) 

Ci = Cfj— {kiT m + Q>)erfc (— — =) /erfc (-^L=) 

+ /«*(£) 


(8.1) 


Here S(t) is the interface position, a is thermal diffusivity, erf is the error function and erfc 
is the complementary error function. This solution presumes a linearized liquidus and solidus 


C s — k s T m , Ci — —kiT m 


(5.2) 


for concentrations at the solid side and the liquid side of the interface, respectively. k s and 
ki are the constant slopes of the straight liquidus and solidus lines. 

There- are two unknown parameters, namely (3 and T m , in the above solution. They 
need to be determined through the heat and mass flux balance conditions at the interface. 
The two interface conditions form a system of two transcendental equations 


pAH/3 = 


K a (T m -Ti)e-£ 2 / a ' | Kl (T m -T 0 
a a y/irerf(l3/a s ) ai y/irerf c(j3/ai) 


(k s 


- ki)T m B = -y[£i- kl - m + Co e -/3 2 /A 
l) mP V 7T erfc (B/VDl) 


(5.3) 

(5.4) 
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It can be proved that there exists a solution for this system. Eqns. (5.3) and (5.4) can be 
solved simultaneously or further reduced to a single transcendental equation 


|a* y/npAH(3a s e rf erfc + %Kia s eii 

(k s - ki)/3e rfc L + ^/Wikie~^ /Dl + 

+ \pDiC^e~^ IDl K s aje - ^ 2/a « + ^ 2 / a * erfc (Jpj 4- Kia s erf (J~^ = 0 


(£)} 


(5.5) 


which determines 0. Then T m is given by 

r wmco£^ 

m [y/n(k s — kl)f3er£c((3/^/D'l) + ^/lTlkle-t 32 / Dl ] 


(5.6) 


In our test we consider a special case in which the following values are chosen for the 
parameters involved, viz. 


p = 1, AH — 1 

K s —Ki = l 

D s = 1, D t = 1/2 

Qig — CLl — 1 

k s = 1 , ki = 2 

T 0 = 1, Ti = -1, Co = 0.1 . 

Substituting (5.7) into (5.5) and (5.6) and solving the nonlinear equations give 
(3 = 0.3526943389 and T m = -0.06872222638 . 


(5.7) 


(5.8) 


Numerical solution is obtained by using FIDAP. Only a finite part of the semi-infinite 
domain, i.e. [0,4], is considered. The boundary conditions imposed for energy equation 
are of mixed type with a constant temperature Ti =— 1 at x=0 and a no-flux (dT/dx = 0) 
condition at x = 4. For concentration, the no-flux condition (dC/dx= 0) is applied to both 
end surfaces. To avoid problems associated with initiation at t = 0, numerical solution 
begins from t = 0.1 with the initial T and C provided by the similarity solution (5.1). 
The backward Euler time integration scheme is used with a time step At = 0.01 and the 
solution ends at t = 0.4. Comparisons are made between the numerical solution and the 
similarity solution. Some typical results are shown in Fig. 5. As we can see the agreement 
is excellent and the error is within 1%. 
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5.2 Axi-Symmetric PoiseuiUe Flow Subject to Buoyancy Force 

This simple test is designed for checking the formulation and FEM implementation on 
the buoyancy force term in the momentum equation. For simplicity, we assume a constant 
change in temperature and concentration which in turn leads to a constant buoyancy body 
force 

p 0 (l - Ar AT - 0 c AC)g . 

Consider a rectangular domain in the axi-symmetric plane with length L and radius 
i?o- The boundary conditions for the momentum equation axe imposed as 

u r — 0, du z /dr = 0 @r = 0 

u r = 0, o z = —A p, @ z = 0 

(5.9) 

u r — 0, u z = 0, @ r = Ro 

u r = 0, a z — 0 . @ z = L 

Furthermore we assume laminar flow and the only non-zero velocity component is u z = u{r). 
For this test problem, the analytical solution can be easily derived and u z has the following 
form 

u(r) = i- a>(1-/3tAT-/3 c AC)s-|| (U§-r 2 ) (5.10) 

where dp/dz is the pressure gradient. The following data are used in our test 

Rq = 0.5, L = 1.0 

PQ = 2, p = 2 

0T = 0.002, p c = - 0.1 (5.11) 

AT = 100, AC = 1 
g = 0.5, dp/dx — —A p = —0.8 . 

The above data lead to a maximum velocity 

“max = “(0) « 0.053125 . (5.12) 

The problem is solved by using the bilinear element. To observe the convergence of 
numerical solution, we start from a coarse mesh with an element size Axi =0.1. We 
then refine the mesh by dividing the square element by half in both coordinate directions. 
Namely, the subsequent meshes have element sizes Ax?, = Ax\ /2 and A ®3 = Axi/4. A 
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Figure 6. The convergence check for FEM 
solution of u max at r = 0 for the axi-symmetric 
Poiseuille flow test problem. The error norm is 
defined as |e| = |u£ ax - u max \, where u£ ax 
is mesh dependent FEM solution and u max is 
the analytical solution. 


verification of the pointwise convergence of u max is presented in Fig. 6. The slope of the 
curve in Fig. 6 indicates a second-order convergence rate for velocity solution, which agrees 
well with the FEM theory [45]. 

5.3 Couette Flow in A Vertical Magnetic Field 

This test is designed for validating the formulation related with magnetic damping. 
Consider the plate-driven simple shearing flow of a Newtonian fluid shown in Fig. 7. The 
flow is subject to a uniform vertical magnetic field. Assume the electric field is everywhere 
equal to zero, the magnetic induction is not distorted by the fluid flow, the magnetic field 
is equal to the imposed value and the buoyancy force term is negligible. 



Figure 7. Plate-driven Couette flow before being subject to a uniform transverse mag- 
netic field. Definition of the problem. 
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In this case, the Lorentz body force is simplified to 

<r m (u x B) x B = — a m u x B 2 i (5.13) 

which acts in the opposite direction to the motion of fluid. Therefore the effect of Lorentz 
force is to reduce the r-component of velocity. Solving the momentum equation with the 
appropriate boundary conditions we obtain the following analytical solution in nondimen- 
sional form 

sinhiHa 

u x = . wrr 

sinh(Ha< 

where Ha is the Hartman number. 

The effect of magnetic damping can be seen more clearly by a very simple asymptotic 
analysis. Consider the limit when Ha becomes very large. It is easy to verify the following 
result 

lim u x = Urn e~ H < h -^U p -4 0 V* < h (5.15) 

ffa-> oo H a— *-oo 

which suggests that when the magnetic field is strong enough the velocity can be diminished 
everywhere in the domain except on the top plate where the above limit approaches U p , 
the moving speed of the plate. Here the effect of the imposed magnetic field is in general 
similar to the effect of increasing the viscosity of the melt. However, care must be taken in 
making analogy between the so-called “magnetic viscosity” with the real viscosity. In some 
aspects, they are quite different. For example, the magnetic viscosity is anisotropic and its 
value depends on the magnitude of the effective velocity component, while the viscosity of 
a Newtonian liquid is isotropic and independent of velocity. 

Due to the continuity of u x , a boundary layer is formed near the top plate and its 
thickness decreases with increase of Ha. This boundary layer is very similar to the so- 
called Hartman boundary layer in the solidifying melt analyzed in [19]. Therefore this 
test problem can also examine how well the thin boundary layer can be resolved in the 
numerical solution. 

In our test the following values are selected for the nondimensional parameters involved 
in (5.14). 

h= 1 , Up- 1, Ha = V2, 10, 50 . (5.16) 

The problem is solved by using the 8 -node 3-D brick element in FIDAP. A total of three 
meshes are used with the element sizes Axi = 1/6, Az 2 = 1/12 and AX 3 = 1/24. A typical 
posterior error analysis based on the three-mesh extrapolation is provided in Table 3. The 




■u Pi 


Uy = U z = 0 


(5.14) 
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Location 

Axi 

Ax 2 

Ax 3 

o 

II 

<1 

Analitical 

x = l/6 

.122772 

.122896 

.122927 

.122936 

.1229367 

II 

h-» 

CO 

.252429 

.252659 

.252716 

.252734 

.2527349 

x = 1/2 

.396241 

.396540 

.396614 

.396638 

.3966391 

x=2/3 

.562272 

.562579 

.562656 

.562680 

.5626810 

x = 5/6 

.759832 

.760054 

.760109 

.760127 

.7601279 


Table 3 Posterior error analysis on numerical solution of velocity by three-mesh extrapolation 
for the magnetic damping Couette flow problem with Ha = \[2. Columns Ax*, (i = 1, 2, 3) 
are the raw u x results obtained from mesh 1 to mesh 3. The values in column Ax = 0 are 
caculated by a quadratic extrapolation. The analytical solution is given by (5.14). 






Figure 8. Numerical solutions for the magnetic damping test problem at higher 
Hartman numbers and comparison with analytical solution. 
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results indicate that the FEM solution does indeed converge at the expected second-order 
convergence rate with mesh refinement and the extrapolated values are accurate to the 
sixth digit. 

To illustrate the difficulties involved with the thin boundary layer in the magnetic 
damping, we solve the test problem at two higher Hartman numbers and present the results 
in Fig. 8. As we can see, at Ha — 10 the coarse mesh is not adequate for providing accurate 
results near the top boundary and one mesh refinement is necessary. When Ha = 50, 
however, mesh 2 with Ax 2 = 1/12 becomes incapable in resolving the thin boundary layer 
and causes the typical numerical oscillation near the top plate. In this case a further 
mesh refinement is necessary in order to obtain more accurate solution. This example 
demonstrates clearly the importance of mesh refinement (or grading) in overcoming the 
numerical difficulties rendered in modeling magnetic damping. 



Figure 9. Numerical solution of the axial solute distribution on the centerline 
and surface of the LTT sample and comparison with the 1-D diffusion controlled 
growth solution. The FEM results are based on the axi-symmetric model for the 
vertical configuration with a constant growth rate V g = lcm/hour. In order to 
display the thin boundary layer, only a portion of the domain is shown. 


5.4 Diffusion-Dominated Growth under 0g 

The last test is a more realistic one and designed for testing the PSSM approach. 
In this test we consider the FEM model given in section 4.1 and compute the numerical 
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solution for an ideal growth condition with no buoyancy-induced convection. To eli minat e 
convection, we set gravity equal zero. It is known from theory that the excess solute 
rejected at the moving interface will be transported by diffusion only. In this case the 
axial solute profile is given by the following 1-D solution [1,2] 

C(z) = Cq ^1 + Vz< 0 (5.17) 

where z = 0 is the interface position. The initial concentration of SnTe in PbSnTe is 
C 0 = 0.2%. The results are presented in Fig. 9. As we can see that the simulated solute 
boundary layer at Op is indeed very close to the diffusion-controlled solution, and hence 
the validity of our numerical model has been verified. 

6. VERTICAL BRIDGMAN GROWTH 

Numerical simulation of Bridgman growth system involves a large number of physical 
and geometric parameters and the cost of each computation has forced us to focus on one 
set of parameters. The particular set of design parameters has been given in Table 1. 
The results presented in this section are therefore based on this set of par am eters and the 
thermal profile given in Fig. 1, unless otherwise stated. 

6.1 Thermal Buoyancy-Induced Convection and Its Effects on Solute Segregation 

In order to study the interaction between thermal and solutal gradients, we first 
consider the density variation caused by thermal expansion only. In this case the thermal 
gradient is the main driving force for the convective flow in the melt. A thermally stable 
configuration (with melt on top) is used in the simulation. Here the varying' parameter 
is the magnitude of gravity ranged from 10“ 4 p to Ip, which leads to a span of thermal 
Rayleigh number from about 7.7 to 7.7 xlO 4 . 

The simulated velocity fields at various gravity levels are plotted in Fig. 10a. Note 
that the plots shown in Fig. 10a are based on the total velocity given in (14) which includes 
a constant translation velocity and the buoyancy-induced convective velocity. At 10~ 4 p 
the velocity vectors shown are straight and uniform flowing toward to the interface, which 
indicates that the buoyancy-driven convection part of the velocity is negligibly small. At 
10 p the convection starts to increase from the center part of the melt and streamlines 
become slightly curved near the interface. However, the convective velocity is still very 
small (about 17% of the growth rate as indicated by the value of « max ). At the gravity 
level of 10 2 g, two flow cells have formed in the melt. The smaller one near the interface is 
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Figure 10a Vector plots of total velocity in the melt for vertical growth at various 
gravity levels. The buoyancy flow is induced by thermal expansion only. The left 
vertical boundary of each plot is the centerline of the sample and the bottom curve 
is the melt-crystal interface. 
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Figure 10b Solute distributions in the melt corresponding to the flow fields shown 
in Fig. 10a. The darker grey scale represents higher solute concentration. 
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driven by the radial temperature gradient. The vector direction in this smaller cell shows 
that the flow adjacent to the interface moves upward at the centerline of the ampoule and 
downward along the ampoule wall. The larger cell is driven by the axial thermal gradient 
caused by the transition in thermal boundary conditions from the adiabatic to hot zone. 
The flow direction in the larger cell is opposite to that in the smaller cell. A third flow cell 
is formed at 10 -l g and the overall strength of buoyancy-induced convection grows with 
the increase of gravity levels. The flow patterns shown in Fig. 10a are in a good qualitative 
agreement with that reported in the literature [19,41]. 

The strength of convection in the melt can be measured by the maximum velocity. 
The effect of different gravity levels on the convective strength in the melt can be seen 
more clearly from Fig. 11 in which the total maximum velocity E/ max versus the magni- 
tude of gravity is plotted for the vertical growth configurations. The convective strength 
level marked by the horizontal dashed line in Fig. 11 is important, because at this level 
the buoyancy-driven convective velocity maximum is the same as the translation velocity 
(growth rate). It has been shown that radial segregation is often maximum when the 
rate of convection and growth are near equal. Therefore it can be considered as a critical 
point for the growth of convective strength in the melt. Let g c denote the gravity level 
corresponding to this critical point. For simulated vertical growth with thermal expansion 
only, g c 5 X 10 4 g. When gravity is less than g c the convection is very weak and grows 
slowly with increases in g. However, when the gravity level exceeds g c convection grows 
concomitantly with g and the effects of gravity on the fluid motion in the melt become 
significant. 

It is interesting to look at the effects of thermal buoyancy-induced convection on the 
solute segregation. There are in general two basic forms of solute segregation, namely 
longitudinal (axial) macro segregation caused by complete mixing and transverse (radial) 
segregation caused by low levels of mixing near the interface. As analyzed by Brown et 
al [19], these two forms of segregation are a function of convection level in the melt. In 
our computation the convective strength is directly related to the magnitude of gravity. 
As shown in Fig. 10b, at low convective flow levels (such as 10~ 4 g and lO -3 #) the radial 
segregation is negligibly small and the axial segregation is very close to the diffusion- 
dominate growth with a thin solute boundary layer near the interface. At 10~ 3 <gr the 
radial segregation starts to increase and the iso-concentration lines become curved. This 
is due to the small convective flow cell formed adjacent to the interface. The excess solute 
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Figure 11. Effects of gravity on the strength of convection in the melt measured 
by the maximum total velocity, U m3uX , which includes the growth rate V g and the 
buoyancy-driven convective velocity. 

rejected at the interface is first swept in a negative radial direction towards the center of 
the ampoule and then is carried away from the interface along the axial direction at the 
center. However, the rejected solute is only carried a very short distance (about twice the 
radius) along the axial direction and cannot go further due to the small size of the flow cell 
driven by the local radial thermal gradient. When the convective flow becomes stronger at 
high gravity levels, the radial segregation continues to increase but the axial segregation 
is still limited within the short range near the interface. Since the rejected solute is only 
circulated in the lower flow cell and cannot be passed to the larger upper flow cell, the 
axial mixing is incomplete. This suggests that for this particular vertical growth condition 
the thermal gradient contributes merely to the low level of mixing near the interface. 

6.2 Effects of Solute Volumetric Expansion 

For non-dilute alloy system, the density variation caused by solute volumetric expan- 
sion becomes important. For PbSnTe, scaling analysis shows that the magnitude of /3 t AT 
and p c AC are globally of the same order. Consequently, the solute buoyancy force is of 
equal importance as the thermal counterpart. In our modeling, the second set of solutions 
considers both the thermal and solute volumetric expansions. The results are presented in 
Figs. 11-13. 
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Figure 12a Contours of stream functions for the vertical growth when both the 
thermal and solutal buoyancy forces are considered. The flow direction of all the flow 
cells, except the top flow cell at 10 ~ 2 g, is clockwise, i.e. the fluid flow upwards along 
the axis and downwards along the wall. Here the calculation of stream function is 
based on the total velocity and 1 ^ = 2.778 \imjsec. 
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Figure 12b Simulated solute segregation in the melt corresponding to the flow fields 
shown in Fig. 12a. The darker grey scale represents higher solute concentration. 



10~ 4 g 10 ~ 3 g 10 ~ 2 g 10 ~ x g 1 g 

Figure 12c The effect of convective flow on the thermal field in the melt for the 
vertical growth when both the thermal and solutal buoyancy forces are considered. 
The values of the isotherm contours are: A = 912.0, B = 958.9, C = 1005.8, 
D = 1052.7 and E= 1099.7° C. 
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The most obvious impact by considering the solute buoyancy force is the changes in 
the convective flow pattern and the significant increase in convective strength. Comparison 
between Fig. 10a and 12a shows that when the solutal expansion effect is included the flow 
field and solute segregation become quite different starting from the gravity level of 10 -3 <?. 
In this case the flow pattern in the melt is greatly changed by the interaction between 
VT and VC. The axial thermal gradient causes material to rise in the center and flow 
down along the wall. This thermal buoyancy induced flow cell (the lower flow cells shown 
in Fig. 10a) is strengthened by the axial concentration gradient which forces the rejected 
lighter solute to be carried away from the interface along the centerline. As a result a 
single large cell is formed at 10~ 3 g. However, the interaction between thermal and solute 
gradient is quite complicated. At 10 ~ 2 g two equal-size flow cells are formed as shown in 
Fig. 12a. The upper flow cell is then diminished again at higher gravity levels. When 
the solute volumetric expansion is included, the strength of the convective flow (measured 
by U m ax) is about one order of magnitude higher than that caused by thermal expansion 
only (for gravity levels greater than g c as shown in Fig. 11). The results presented here 
demonstrate the significance of the effect of solute volumetric expansion and the interaction 
between the solutal and thermal fields for non-dilute alloys. 

It is interesting to observe how the convective flow field changes with increase of 
gravity levels. At low gravity level, such as lO -3 #, the flow cell shown in Fig. 12a indicates 
that the flow at the lower part of the ampoule can be divided vertically into two opposite 
flow regions, with the division fine lies approximately at the middle of the radius, i.e. 
r = Q.5R. The central part r < 0.5 R is a cylinder with an upward flow, while the outer 
annular region 0.5i? <r <R has a downward flow. Since the cross-section areas of these 
two regions are unequal, the average flow speed on them is different. In order to satisfy 
the continuity condition, the flow is faster in the smaller central region than that in the 
larger outer annular region. However, with increase of convection strength, the flow cell 
is pushed towards to the wall. As we can see from Fig. 12a, the outer annu lar region 
becomes smaller and the division line at lg is much closer to the wall. This results in a 
thin boundary layer and a faster flow near the wall at high gravity levels. The thermal 
contour plots in Fig. 12c also indicate a thin thermal boundary layer at lg level. These 
thin boundary layers render great difficulties for numerical solution, causing slow-down in 
convergence or even divergence of the nonlinear iterations. These difficulties occurred in 
our computation were also reported in the literature [19]. 
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The effect of solute gradient on segregation is also significant. On the interface, the 
radial solute gradient acts in an opposite direction to the thermal radial gradient and tends 
to reduce the radial segregation. At the center of the sample, however, the axial solute 
gradient acts in the same direction as the axial thermal gradient and tends to enhance the 
axial segregation. This results in a weaker local radial segregation near the interface and 
a much more complete mixing at high gravity levels, as can be seen in Fig. 12b and 13. 

In order to describe quantitatively the transition from diffusion-controlled growth 
(without bulk convection) to growth with intensive laminar convective mixing, we compute 
two conventional segregation parameters. The first is the effective segregation coefficient k e fi 
defined as [19] 

fc eff = fc<C>j/<(7» (6.1) 

where < C >/ is the lateral average of solute concentration over the solid/liquid interface 
and is the volumetric average of solute concentration in the melt. The second is 

the percentage radial segregation AC defined as 

AC = \5C\ max / Co (6.2) 

or 

AC = |C t0 p — ^bottom |/ Cq . (6-3) 

Here |^C| max is the maximum difference in concentration at the interface. For directional 
solidification, diffusion-controlled growth with a planar interface leads to uniform solute 
distribution in radial direction, thus we have AC = 0. If the melt is sufficiently long that 
the diffusion layer near the interface occupies only a small fraction of the total length, fe e ff 
approaches unity when there is no convective flow in the melt other than the unidirectional 
growth velocity V^. The other limit of fc e f? represents the steady state well-mixed growth 
in which the intense convection leads to a complete mixing in the melt. In this case the 
value of fc e ff approaches k. 

The computed segregation coefficients for the vertical growth are plotted as a function 
of gravity levels in Fig. 13. Here the quantitative axial and radial segregation analysis is 
consistent with the results seen in Fig. 10b and 12b. For thermal buoyancy induced flow 
case, the values of fc e ff are very close to one at all gravity levels, which indicates the 
diffusion-dominated growth in axial segregation even though the local radial segregation 
may be quite large at high gravity levels as suggested by the values of AC. 
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k = 0.7 


Figure 13. A quantitative description of the effects of gravity on the axial solute segrega- 
tion (a) measured in terms of the effective segregation coefficient, fc e ff and the radial solute 
segregation (b). The results axe based on the axisymmetric model for the vertical growth and 
the values of AC were calculated from the liquid side of the interface. The results labeled 
by fit considers thermal buoyancy force only. 
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When both thermal and solute buoyancy forces are considered, the calculated fc e ff 
values in Fig. 13(a) show a diffusion-controlled growth at low gravity levels. k eS is very 
close to k at high gravity levels which indicates that a well-mixed state has been achieved. 
Radial segregation results show that weak convection (in the range 10 — 6 g to 10~ 4 g) and 
strong convection (at lg) both produce low radial segregation, while the intermediate 
gravity levels produce high radial segregation. As shown in Fig. 13(b), AC has a local 
mavimiim at 1O _3 0, which is very close to the gravity level predicted by the critical g c 
shown in Fig. 11. The plots of fc eff and AC v.s. the magnitude of gravity presented here 
agree well with the schematic curves given by Kim, Adomato and Brown in [19]. 



B (KG) 


Figure 14. Effects of an axially imposed magnetic field on the strength of 
convective flow in the melt measured by the maxim u m total velocity, U max . The 
results are based on the vertical growth model with R = 1 cm, = 0.5 pm j sec 

and Pc = 0. 


6.3 Magnetic Damping 

The action of an axial magnetic field on convection in the melt is to interfere with 
the radial velocity component. As pointed out in ref. [19], at large B the magnitude 
of the radial velocity decreases proportionally to the square of the strength of magnetic 
induction. The axial velocity component is not affected by the magnetic field except in 
the coupling through the incompressibility condition. When the magnetic field is strong 
enough, an almost uniaxial flow (i.e. u z = V g , u r = 0) can always be obtained. One of the 
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goals of this study is to quantitatively determine the required magnetic field levels at which 
the diffusion-controlled growth can be achieved in orbit. In Fig. 14 and 15 we provide a 
typical analysis of the effects of magnetic damping for vertical Bridgman at a low gravity 
of 10“V Note that the results presented in this subsection are based on a larger ampoule 
radius R = 1cm and a slower growth rate V g = 0.5 fim/sec. With these conditions, 10 ~ 5 g 
resulted in an almost diffusive growth without magnetic field. 

Our numerical results suggest that magnetic damping on the convective flow in the 
melt is practically effective for micro-gravity crystal growth. By considering the magnetic 
damping as a continuous loading process with the increase of B, the typical Z7 max v.s. B 
curve given in Fig. 14 shows that there are mainly two stages during the magnetic damping 
process. The first is a rapid change stage in which the strength of the convective flow in 
the melt is significantly reduced with an increase of magnetic field strength. The second is 
a slow (asymptotic) change stage during which the velocity decreases very slowly with an 
increase of B. These two stages can be clearly seen in Fig. 14. In this case, over 90% of the 
original strength of convection is damped when B is increased from zero to 3000 Gauss. 
However it requires another 6000 Gauss to eliminate nearly all of the original velocity 
strength. 

In Fig. 15 we plot the interface shape, radial and axial solute segregation at several 
B levels. Fig. 15(b) shows how magnetic damping affects the radial segregation on the 
interface. Here the change of radial segregation with the increase of B is not monotonic, 
it is reduced to a local minimum at about B = 1.5KG and then increased again. Fig. 15(c) 
and (d) shows the effectiveness of magnetic damping on the axial segregation under the 
microgravity condition. At B « ZKG the axial segregation is very close to a diffusion- 
dominated growth state. 

At high gravity levels, our results suggest that the axially imposed magnetic field is 
not effective for suppressing convective flow in the melt. Fig. 16 shows an example of 
magnetic damping for the vertical growth on earth. As we can see that the u max is still 
about 1000 times higher than V g at B = 25 AG and the decrease of u max is much slower 
than that at the gravity level of 10~ 2 ^. The scaling analysis suggests that the req uir ed 
magnetic field strength for the terrestrial growth will be much higher than the generally 
available 50 KG in order to achieve the diffusion-dominated growth condition. 
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B =0 


B = 1.538KG 
B = 3.077KG 
■ B = 8.769KG 


B = 0 

B = 1.538KG 
B = 3.077KG 
No Convection (Og) 


B =0 

B = 1 .538KG 
B =3.07710? 

No Convection (Og) 


Figure 15. Numerical solutions of the vertical Bridgman growth with increasing magnetic 
field strength measured in kilo-Gauss. (a) Liquid/solid interface shapes; (b) Radial solute 
segregation at interface; (c) Axial segregation at the center of ampoule (R = 0); (d) Axial 
segregation on the surface of sample (R= 1.0cm). 
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Figure 16. Magnetic damping for the high gravity level of 1 g. The results 
are based on the vertical growth model with R = 1cm, V„ — 0. Sum/ sec and 
&= 0 . 
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7. HORIZONTAL BRIDGMAN GROWTH 

Vertical Bridgman represents an idealized configuration which may only exist on earth. 
In space, the situation is quite different. First, it is practically impossible to exactly align 
the furnace axis with the residual gravitational vector. Second, any deviation of the furnace 
axis from gravity direction will dramatically change the flow field. The results presented 
in Fig. 17(a) and (b) suggest that even a 1° mis-alignment at low gravity of 10 ~ 4 g can 
completely distort the axisymmetric flow pattern and increase the strength of convection 
by a factor of up to seven in our particular case. 

==1 ' ~ ^ ~ ' — =§ 10 ~ A g 

— — — — ► 

(a) Axi-Symmetric Solution of Vertical Bridgman U m&x = 1.17V* 

1 ° 

10 -4 5 

(b) 2-D Solution with Mis-Aligned Gravity U m&x — 8.25V S 

1 

10 ~ 4 g { 



I . 

10~ 2 5 , 

(e) Horizontal Bridgman U m ax = 1061.3V* 

Figure 17. Effects of orientation and magnitude of gravity vector on the fluid motion in 
the melt of PbSnTe. 





To study the effects of the orientation and magnitude of gravity, we present in 
Fig. 17(c)-(d) the flow fields in the melt, and in Fig. 18 the interface shape, radial and 
axial segregation for horizontal Bridgman growth at four typical gravity levels. Note That 
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the solutions presented in this section consider the density variation caused by both the 
thermal and solute volumetric expansions. 

For horizontal Bridgman, the growth at 10 -5 <? is very close to diffusion-dominated 
state as indicated in the axial segregation plots in Fig. 18(c) and (d). At 10~ 4 <? the thermal 
buoyancy induced flow cell in the center part rotates in the same direction as in the vertical 
Bridgman configuration. However, the radial solute segregation and VC y are locally strong 
enough to force the rejected solute flow up along the interface and form a clockwise flow 
cell near the interface. For a given growth speed V g , this solutal buoyancy induced flow 
cell exists only at low g when the mixing is incomplete and becomes less pronounced as 
gravity increases. When gravity > 10~ 2 ^ the small flow cell near the interface vanishes 
and the thermal convection grows to a larger cell. This strong convection results in a 
complete mixing and a much more uniform solute distribution along both the axial and 
radial directions, as can be seen in Fig. 18(b)-(d). 

For the horizontal growth, the axial magnetic field is less efficient than for the vertical 
growth. A typical magnetic damping analysis is presented in Fig. 19. As we can see from 
Fig. 19(d), a three Kilo-Gauss field can bring down about 90% of the convection strength 
at 10~V However, a much higher magnetic field (about 40 KG) is required to reach to 
the diffusion-dominated growth state. 

8. CONCLUSIONS 

Since the alignment of furnace axis with residual gravitational acceleration aboard any 
orbiter cannot be practically guaranteed to within a small fraction of a degree, we should 
only resort to examining the results obtained from the horizontal Bridgman growth models. 
Even though the DC acceleration levels typical on the shuttle are < 10 ~ 5 g (depending on 
the furnace location relative to C.G.), we recommend using a worst case scenario in which 
the ^-vector is perpendicular to the ampoule axis. Furthermore, it has been known that 
isolated large amplitude peaks in the order of O(10 -3 )p are common in space. Therefore, 
a conservative estimate at this time can probably be obtained by using the horizontal 
Bridgman model with 10~ 4 p DC case. Our results show that a magnetic field of few kilo- 
Gauss will be sufficient to substantially reduce the effects of buoyancy-induced convection 
on the solute segregation. Investigation of a 10 ~ 5 g DC case with a time-dependent 10 _3 p 
acceleration pulse will be pursued later. 

Future plans call for ground work in a magnetic furnace and continued modeling 
work on predicting and comparing with experimental data for a systematic examination 
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Figure 18. Modeled horizontal Bridgman growth at V g = lcm/hour and three low gravity 
levels, (a) Liquid/solid interface shapes; (b) Radial solute segregation at interface; (c) Axial 
segregation on bottom boundary ( y = —0.5cm); (d) Axial segregation on top (y = 0.5cm). 
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Figure 19. Numerical solutions of the horizontal Bridgman growth with increasing magnetic 
field strength measured in kilo-Gauss. (a) Liquid/solid interface shapes; (b) Radial solute 
segregation at interface; (c) Axial segregation on bottom boundary ( y = —0.5cm); (d) 
Damping of U max with increase of B. 
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of the phenomena at hand, and for specifying the required magnetic field for microgravity 
experiments. 
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